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Derivation of the Ewald Summation Formula. 


The total Coulomb energy of a system of N particles in a cubic box of size L and their 
infinite replicas in PBC is given as: 
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where q; is the charge of particle i. The main cell is located at n = (0,0, 0). Its image cells 
are located at L,, intervals in all three dimensions as n goes to infinity. 


The unit cell vector is n = (nj,n2,n3) = 1yLk + nzLV + n3LZ where %, ¥, Z are the 
Cartesian coordinate unit vectors. The first sum in Coulomb energy equation is primed to 
indicate that terms with i = j are omitted when n = 0. 


This is a single and conditionally convergent series which means that the result depends 
on the order of summation. Here the q; and q; represents the two-point charges with 
some constant values. From the nature of this representation we will observe that the 
amount of consumed time will increase geometrically with the increase in the number of 
particles in the system. To decrease the computation time, the charge distribution notion 
is used instead of point charge notion. This had enabled Ewald and Bertaut to develop an 
absolutely convergent series to compute Coulombic interactions. 


Bertaut, replacing the ionic crystal by a spherically symmetric charge distribution, 
calculated the Coulombic interaction energy by a single absolutely convergent infinite 
series whose terms are a function of the lattice vectors in the reciprocal space. Bertaut's 
series is absolutely convergent but is not very efficient if computation time is taken as the 
parameter. 


The disadvantage of Coulomb's energy equation is that basically it is a single and 
conditionally convergent series. These two characteristics of this formulation form an 
incompatibility for high speed computing environments. 


Ewald Summation provided an efficient technique to evaluate the long-range 
interactions between particles and their periodic images. Ewald Summation re- 
formulates the Coulomb's equation with the sum of two rapidly converging series plus a 
constant term. These series calculate the interactions at different ranges and the 
summation of these terms gives us the Coulomb energy within an error range. The Ewald 
Summation is formulated as: 
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where @, is the direct (real) sum, ¢,. is the reciprocal (imaginary) sum and @¢, is the 
constant term. The problem is to calculate the electrostatic potential experienced by one 
ion in the presence of all the other ions in the crystal. We consider a lattice made up of 
ions with positive or negative charges and shall assume that the ions are spherical. We 
compute the total potential ¢ at an ion as the sum of two distinct but related potentials 
(excluding the constant term). 


Force expressions can be easily obtained from Ewald energy formulation by direct 
differentiation in each coordinate p = x,y,z. To avoid complexity, the equations in this 
paper will be potential energy based. And the self-term will not be shown since it, as a 
constant, will disappear in force calculations. The equations for these terms for the 
Gaussian distribution are, for the direct sum: 
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where €, is the permitivity of free space, a is the Ewald coefficient, N is the number of 
charges in the system and 
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The reciprocal sum: 
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where V is the volume of the simulation box, m is the reciprocal space vector. 


In the reciprocal sum part of the E.S., each point charge is replaced by a charge 
distribution of proportional magnitude and equal sign which spreads out from the lattice 
site in a spherically symmetrical way. Or in other words each ionic charge is distributed 
within a sphere of radius R so that the charge density is proportional to a normalized 
distribution function y(r), where r is the distance from the atomic center. In the direct 
sum part each particle is represented by the original point charge and a charge 
distribution of proportional magnitude and opposite sign. We see the significance of the 
shape of the charge distribution on the computation time clearly at Ewald Summation. 
By changing the width of the Gaussian type charge distribution (or the Ewald coefficient 
parameter) we are able the shift the amount of work done between the direct sum and 
the reciprocal sum parts. As the distribution gets narrower, the reciprocal sum becomes 
more dominant in the result. At the far end, if the distribution function is selected to be a 
Dirac function then the point charge and the distribution function cancels out at the 
direct sum and the result is generated only by the reciprocal sum. In general, this work 
load shifting aims to prevent a bottleneck that may come from the computation of one 
term; either direct or reciprocal sum parts. A lot of work has been done to find the 
optimum Ewald coefficient parameter. 


It is worth pointing out that in distributed environments determining the optimum 
Ewald coefficient parameter should consider the structure of the distributed 
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environment too; which adds the problem of finding an optimal value a new perspective. 
The size of the processing system, the structure of processing elements, the distribution 
algorithms, workload allocation techniques implemented are some factors that affect the 
optimal value. 


An electrostatic potential is a long-range function, which decays only as the inverse 


power of the distance. An alternative to the direct application is to split the potential into 
three parts as mentioned before: 
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, where (r) is the total potential of the system, ,(r) is the direct sum potential, ,(r) is 
the reciprocal sum potential and ,(r) is the correction term. 


The total potential can be represented in terms of electric field E : 
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where y(r) represents the shape of the charge distribution as a function of distance, r. 
The evaluation of this term for Gaussian shape charge distribution below: 
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If the charge distribution is zero at some distance d,(or considered to be zero) then the 
electrostatic potential of the y(r) shaped charge will be equal to the electrostatic potential 
of a point charge qx for r > d. Since the ¢,(r) is the difference of these potentials, it can 
be written as: 


dba(r) = O(r) — $- (7) 
For Gaussian distribution: 
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In the case of a finite distribution (e.g. triangular shape charge distribution) the upper 
limit of the second integral in equation should be the point where the distribution hits 
zero (e.g. 1/a). 


The derivation of the reciprocal sum part for a particular charge distribution involves 
computing the Fourier transformation of the charge distribution. 


The Fourier transform of a charge distribution, y(r) is: 
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which for Gaussian distribution yields: 
¥(k) = exp(—k*11*/a) 


With the reciprocal sum potential 
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For Gaussian distribution we get: 
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